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Abstract 



We consider a d-dimensional crystal with an arbitrary harmonic interaction and an an- 
harmonic on-site potential, with stochastic Langevin heat bath at each site. We develop an 
j^ I integral formalism for the correlation functions that is suitable for the study of their relax- 

Ch ' ation (time decay) as well as their behavior in space. Furthermore, in a perturbative analysis, 

I I for the one-dimensional system with weak coupling between the sites and small quartic an- 

' O ' harmonicity, we investigate the steady state and show that the Fourier's law holds. We also 

obtain an expression for the thermal conductivity (for arbitrary next-neighbor interactions) 
and give the temperature profile in the steady state. 



PACS: 05.70.Ln; 05.40.-a; 05.45.-a; 44.10. -fi 
'^ I Short title: Anharmonic Crystals with Self-consistent Reservoirs 

^ ! 1 Introduction 

^: 

^D ■ We are surrounded by phenomena involving non-equilibrium processes, but our understanding 

of such systems, i.e., the number of models that permit detailed calculations, is very limited. 
In particular, it is unknown a simple way of finding the properties in the steady states: e.g., 
a rigorous derivation of the (phenomenological) Fourier's law from a microscopic anharmonic 
'^ i Hamiltonian model has not been established up to now (see pQ, [2] for a review). It makes the 

Q I analysis of simple dynamical models describing non-equilibrium processes a problem of interest. 

O ' A commonly studied microscopic model is the Hamiltonian chain (or its d-dimensional ver- 

sion) of N interacting oscillators coupled to heat baths at each site or at the boundaries only, 
K^ \ and its anharmonic version with small quartic on-site interactions. 

^ • For the harmonic case of the model with thermal reservoirs at the boundaries, the covariance 



of the stationary state has been calculated in [S|, long time ago. There, it is shown that the 
heat current is independent on the length of the chain, and so, the Fourier's law does not hold. 
The temperature profile is also computed in [3|: the temperature is essencially constant in the 
interior of the chain, but decreases exponentially close to the hotter bath and increases close to 
the opposite end. I.e., the profile has the lowest temperature near the hottest reservoir and the 
highest temperature near the coldest reservoir. For the anharmonic case, there are interesting 
and recent results. The existence of steady states is proved in ^, and the positivity of entropy 



production in [Sj. Numerical results strongly suggest that the Fourier's law holds in such a case 
jH], HI, but, in the opposite direction, a perturbative analysis jHI shows that the heat current 
does not depend on the size of the system. Also in the perturbative study 8^, as in harmonic 
case jS], the temperature profile (discarding the exponential decay in the bulk of the chain) is 
in the "wrong" way: the hottest temperature is near the coldest bath, and vice-versa. In short, 
it is unclear whether the Fourier's law holds or not in such anharmonic models. It is worth to 
recall that other results also indicate that it is wrong the opinion that the sole anharmonicity of 
the on-site potential shall ensure normal heat conductivity in some commonly used models jSj- 

The harmonic crystal model with next-neighbor interactions and heat bath at each site has 
been recently analyzed in ^01- It is proved, for a uniquely fixed temperature profile leading 
to the steady state (given the temperatures at the boundaries), that the heat current satisfies 
the Fourier's law. For the case of more intricate interactions (intense and beyond next-neighbor 
sites), for a chain with some few sites, some results presented in jllj indicate that there is a 
"strange" heat flux in the harmonic network (and the authors claim that the results persist 
under weak anharmonic perturbations): inside the chain, the direction of the heat fluxes can 
not (in general) be supposed from the temperature of the heat baths. 

In the present paper, also with the aim of studying the dynamics of simple microscopic models 
in order to understand properties of non-equilibrium systems, we study the anharmonic version 
of this crystal with stochastic Langevin heat bath at each site (model named as crystal with self- 
consistent reservoirs). We describe an approach and obtain an integral formalism suitable for the 
study of the correlation functions (of the d-dimensional system with quite general interactions). 
Furthermore, using perturbative calculations, for a weak coupling between the sites and a weak 
anharmonic potential, we show (for the one-dimensional system) that the Fourier's law still holds. 
That is, we show (at least up to first order in the perturbative computation) that the Fourier's 
law is valid for this microscopic anharmonic Hamiltonian model. We also obtain an expression 
for the thermal conductivity (for next-neighbor interactions which may arbitrarily change along 
the chain), and give the temperature profile in the steady state. For the simpler case of next- 
neighbor interactions constant along the chain, our results (considering the anharmonic model) 
coincide with those of the harmonic case recently described in [IDj . 

The rest of the paper is organized as follows. In section 2 we present the model and some 
expressions for the energy current. The integral formalism for the correlation functions is devel- 
oped in section 3. In section 4, in a perturbative computation, we analyze the energy current in 
the steady state and the Fourier's law. In section 5 we argue on the reliability of the perturbative 
results and present some concluding remarks. 

2 The Model and Initial Considerations 

Lets us introduce the model to be analyzed here and some expressions for the energy current. 
We consider the stochastic Langevin dynamics of an anharmonic crystal, i.e. a scalar field 
lattice model with unbounded spin variables in a d-dimensional lattice space box A C Z , with 
stochastic heat bath at each site. Precisely, we take a system of A^ oscillators with Hamiltonian 

N ^ N N 

j=i j^i=i j=i 



where M > , V gives the anharmonic on-site perturbation (e.g., 'P{qj) = Qj), and we consider 
the time evolution given by the stochastic differential equations 

dqj = Pjdt, j = l,...,iV, (2) 

dpj = -g—dt - CPjdt + -ty^dBj, j = 1,...,N, 

where Bj are independent Wiener processes, i.e., dBj/dt are independent white noises; C is the 
heat bath coupling; and jj = 2(^Tj, where Tj is the temperature of the j-th heat bath. 

To describe the energy current in the system, we write the local energy of the spin (oscillator) 
j as 

H,iq,p) = \p] + U^'\q,) + IY1 U'^'^il^ - '?')' (3) 

where the expression for U^^' and U^'^' follow immediately from (^ and X^,=x ^j ~ ^- Then, 
we have 



dHj{t) 



{Rj{t))-{r,<,-r,y), (4) 



dt 
where (•) denotes the expectation with respect to the noise distribution, and 

{R,{t)) = Q {T, - {p])) (5) 

gives the energy flux from the j-th reservoir to the j-th site. The energy current inside the system 
is given by J-j, where 

In particular, in the steady state we have {dHj{t)/dt) = 0. We will turn to these expressions to 
discuss the Fourier's law later. 

3 The Integral Formalism for the Correlation Functions 

For convenience, we introduce the phase-space vector (j) = (g, p) with 2N coordinates and write 
the equation for the dynamics (2) as 

4) = -Acp - XP'{<i)) + ar], (7) 

where A = [AP + J) and a are 2N x 2N matrices given by 

I above is the unit N x N matrix; J is the N x N matrix for the two site interaction Jij 
(see(pP)); and Ai,T,T are diagonal N x N matrices: Aiji = M6ji , Tji = (6ji , Tji = Tj6ji. rj are 
independent white-noises; V'{(f)) is a 2N x 1 matrix with V'{(j))j = for j = 1, . . . ,N and 

p'(^), = ^^^izi^ for i = N + l,...,2N. (9) 

d(Pi-N 



To describe the dynamics we first consider the system without the couphng J among the 
sites and without the anharmonic perturbation (A = 0) (interactions which we include in a 
second step). Then the (straightforward) solution of above with J = 0, A = is the 
Ornstein-Uhlenbeck process given by 



^(t) = e-*^V(0) + / rfse-(*-^)^V7?(s). (10) 

Jo 

For simplicity we take (j){0) = 0. The covariance of this Gaussian process evolves as 




sAO-sAo'^ 



C{t,t)= / dse-^^ Ve-^^ . (12) 

Jo 

It is easy to see (e.g. diagonalizing ^4*^ ) that 

expH^«)=.-ic„shM{/ + i?=y^( j[^ _-^^)}. (13) 

(and a similar expression follows for the transpose exp{—tA^ ) ), where I is the N x N unit 
matrix, etc; p = ((C/2)^ — M) (we assume that (C/2) > M > 0). In this simple case (of 
J = 0, A = 0), as t ^ oo we have a convergence to equilibrium (any single site is isolate) and 
the stationary state is Gaussian, with mean zero and covariance 

C = fds e-^^V-^^°" = ( I °), (14) 

where, again, T is a diagonal matrix with elements Ti6ij (in short, for any site we have a Gibbs 
measure at temperature Tj). 

To introduce the coupling interactions and the anharmonic potential, we use a tool of general 
theory of stochastic differential equations, namely the Girsanov theorem ^2- It gives a measure 
p for the new process (jT)) as a "perturbation" of the measure pc associated to the process with 
J = , A = 0. Precisely, for any measurable set A, it states that p{A) = Eo{lAZ{t)), where 
Eq is the expectation for pc (the process with J = , A = 0); 1a denotes the characteristic 
function, and 

Z{t) = exp( u-dB-- u^dsj, (15) 

7i Ui = -Jik(t)k-yP'{(j))i, 

(the inner products above are in M^^). From Q and the expression above for Ui, note that Ui 
is nonvanishing only for i > N (i.e., i £ [A^ + 1, A^ + 2, . . . , 2N]). In what follows we will use the 
following index notation: i for index values in the set [A^ + 1, A^ + 2, . . . , 2A^], j for values in the 
set [1,2,...,N], and k for values in [1,2,..., 2A^]. 

For clearness, let us rewrite the stochastic equations for the initial process (with J = , 
A = 0) as 

dcPj = -A%cf)kdt, j £[!,..., N], (16) 

#i = -A^i^cPkdt + %^^dB,, i€[N + l,...,2N], 



the sum over k (in [1,2..., 2A^]) is assumed above (as well as obvious sum over some indices in 
what follows). 

Turning to the terms in Z{t) we have 

UidBi = 7^ Uj • 7- ' dBi = 7^" Uj • [d4>i + A'^,^4>kdt) 

that follows from H15|) and (jlbf) above. We still use Ito formula to write the terms with dcpi as 

--il^Jij(t)jd4)i = -dFi - -^^^4)iJijA^jf.(t)kdt, 

With similar manipulations we obtain 



Z{t) = exp I u ■ dB — - u ds 



exp{-Fi(0(t)) + Fi(0(O)) - AF2(,/.(t)) + AF2(</.(0))} 

X exp |- / Wj{4>{s))ds - [ XWx{cP{s))ds - I XWxj{(t){s))ds \ , (17) 



./O JO 



with 



Fiim) = l~^4>i{t)JiMi), F2{4>{t)) = li^V\(t>)i{t)(t>i{t) 

+-</>i'(s)j/i7r'^ij<^j(^)> 
\Wx{<P{s)) = \-i-'Us)V"{<i^Us)AtN,k<t^k{s) + \-ii'V'{<PUs)A%<l,k{s) + 

+\\h;\V'{<l^).f{s)At^^,<t>k{s). 

And so, for the expectations, considering the process with coupling between sites and anharmonic 
perturbation, we have, e.g. for the two-point function, 

{Mh)(t>v{t2)) = f Mh)Mt2)Z{t)dfic{(t>), tiM < t. (18) 

The formula above, a Feynman-Kac type integral representation, is suitable for the study of 
general n-point correlation functions: for the analysis of their time decay (relaxation properties), 
space behavior, etc. In particular, we will analyze the energy current in the steady state, problem 
that involves the investigation of terms such as limt^oo {4'iii)4'ji't)): see (jH)) . 

4 The Heat Flow and Fourier's Law 

To study the heat flow in the steady state we need to analyze the two-point correlation functions 
given by formula ©. The averages over the stationary distributions will be obtained as the limit 

{4>u(t>v) = lim {(l)u{t)4>v{t)) = lim / 4)u{t)(l)v{t)Z{t)diic{(t>)- 

t— >oo I— >oo / 




We will establish conditions for the convergence to the steady state later. 

To carry out the computation, note that C{t,s), given by ((TT |) - ((TH) . may be written as (for 
t> s) 

C{t, s) = exp(-(t - s)A'^)C + O (exp[-(t + s)C/2]) , 

and the effects of the second term in the r.h.s of the equation above disappears in the correlation 
formula in the limit of t ^ oo. 

For the anharmonic interaction we choose, for ease of computation, V{4>)i{s) = ^ : (pf_j^{s) :, 
where the dots mean the Wick order with respect to the Gaussian measure fie- 

We will make a perturbative analysis , i.e., we will assume that the coupling between two 
sites J as well as the anharmonic potential coefficient A are small. Hence, up to first order in J 
and A, after (considerable but straightforward) calculations we have 



2^7 [Jv+n,u-nTu~.n - Ju,vTv] 5u-N,v for u £ [N + 1, . . . , 2N],v £ [1, . . . , iV] , 

(19) 
Tu5u,v for u,v£ [iV + l,...,2iV]. 

For simplicity we will restrict the analysis of the energy current to one-dimensional systems 
only. From © we have 

^,<= E J^^NA^.-^r)^-^^^^^^^^^, (20) 

1" > j 

re[l,...,N] 

where {J'j<) denotes de energy flow between site j and the site r (with r > j) connected by the 
interaction J . Using the results describe in (|T9|) above we obtain 

^^k) = E%^(^-^.). (21) 

r>j 

Let us analyze, in particular, the case of next-neighbor interaction only. In such a case 

:f,^,^, ^ CF,^) = ^^g^ (r,+i - r,) . (22) 

The condition (dHi/dt) = 0, that characterizes the stationary state , together with expres- 
sions (|1I5)) and {Rj{t)) = (that comes from © and (UH)) lead to 

•^1^2 = -^2^3 = -^3^4 = . . . = JFjv-i^AT. (23) 

I.e., using the notation Jj = (J'j^jyj^i) /2(M, 

Ji (T2 - Ti) = J2 (Ts - Ta) = Jg (r4 - Tg) = . . . = Jn-i {Tn - T^-i) • (24) 

It is easy to see that given the temperatures at the boundaries Ti and T^, and nonvanish- 
ing Ji, J2, . . . , Jn-1, there exists an unique solution T2,T3, . . . ,T/v-i for the linear system of 
equations above ((23}. Namely, we obtain 

11 1 \"^ / 1 1 1 



^■ + (x + ^ + - + *:7 nX^X^-^T^K^"-^''- P5) 



that determines the temperature profile in the steady state. Note that it is a monotonia function, 
oriented in the "right" way: the hottest temperature is near the hottest bath, and vice-versa. 
For the energy current we get 



A^-1 ' A^-1 \Ji J2 Jn-1 

' (26) 

that is, the Fourier's law still holds. For the simpler case of the same interaction between two 
any next-neighbor sites, i.e. Ji = J2 = . . . = Jn-i, we have, for the thermal conductivity, 

X = Ji = ^^^^- (27) 

For comparison, in jlUj the authors treat the linear dynamical problem, i.e. Q with A = 
and 

V^ C J' =LO^-6r+l,j-6r^ij + {2 + U^)6r,j), ^^^ 

and obtain (in a nonperturbative approach) 

^2 1 
y = — f29) 

In our case (considering the same J of |lUj). we have 

A-( ^ -^ \ "^ =UJ^i-5r+l,j-5r-l,j), , . 



J + M C y' M =M6rj, M = uj\2 + u 

and so, our formula (j^TJ above becomes 

{-u^'f _ ^' .... 

^ 2C(2 + z^2)^2 ((4 + 21/2)- y-^ ) 

Considering that our computation was carried out in a perturbative approach with small J but 
M not small (see the comments at the final section ), i.e., uj'^ small, u large, we have in (29) 



aM(4T^«^'(i + ^^) 



1/2 + 2. 



That is, our computation, when restricted to the case treated in jTUI, leads to the same result. 
In short, we have shown (in a perturbative analysis: up to first order in A and J) that the 
Fourier's law still holds for the harmonic crystal with self-consistent reservoirs when a small 
nonharmonic on site perturbation is introduced in the interaction. 

5 Concluding Remarks 

The approach presented here establishes an integral representation for the correlation functions, 
say, a Feynman-Kac type formalism. That is, in some sense, we map the stochastic problem on 
a non-canonical field theory. Such an approach is inspired on previous works considering the 



study of the relaxation to equilibrium of some nonconservative stochastic Langevin systems ^31 ; 
[Tl] . jT3], jTB], |T7]. There, the time decay of the two and four-point functions is analyzed in 
detail. A perturbative study is carried out within the integral formalism in the regimes of low 
and high temperature. In the low temperature region, for the system with a weak anharmonic 
potential and a bare mass (the coefficient of the local quadratic term) large enough, it is proved 
that the perturbative analysis is not naive: e.g., the rigorous results described in ^3] show that 
the complete treatment of the four-point function adds only small corrections to the behavior 
obtained by the perturbative calculations presented in ^Hj. Using similar techniques (cluster 
expansions, etc) we expect to prove the results about the behavior of the correlations presented 
here (for small A, nonzero M and C,, and Tj not large). The perturbative analysis of our system 
with all the reservoirs at (different but) high temperature (i.e., with the perturbative parameter 
given by 1/Tj instead of A) shall be possible following procedures similar to those described in 
|16j and references there in. 

Another interesting open problem is the behavior of the system in the limit of the coupling 
with the interior heat bath taken to zero: in such a case, as we have mentioned before (compare 
and with [Hj), it is not clear if the Fourier's law is valid or not. 
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